function ray_make_traveltime(varargin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function ray_make_traveltime
%
% This function is internal to the ray_invert_db and ray_invert_ps functions.
% It generates the traveltimes.tsv file required by raytrace3d invert and 
% raytrace3d PSinvert. For S-waves, the function assumes a constant P:S 
% ratio of 1.73.
%
% This function requires the presence of one other function:
%   ray_latlon2xyz.m
%
% Required Inputs:
%   db:                 A database pointer containing, at a minimum, the 
%                       site, arrival, assoc, and origin tables.
%
%
% Optional Input:
%   locations:          The path to the output file containing updated
%                       earthquake hypocenters generated by the ray_locate 
%                       function
%
% Output (to filesystem):
%   traveltimes.tsv:    A file containing the station-origin information
%                       required by the raytrace3d invert routine
%
% Author: 
% Matt Gardine
% June 2009
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if exist('ray_defaults','file')==2
    [ref_lat,ref_lon,ref_alt,projection]=ray_defaults();
else
    ref_lat=17.01;
    ref_lon=-105.99;
    ref_alt=0;
    projection='flat';
end

fid = fopen('./traveltimes.tsv','wt');

switch nargin
    case 1
        db = varargin{1};
        
        db = dbsort(db,'orid');
        
        [orig_lat,orig_lon,orig_depth,sta_lat,sta_lon,sta_elev,orig_time,arr_time,phase]=...
            dbgetv(db,'origin.lat','origin.lon','origin.depth','site.lat','site.lon','site.elev','origin.time','arrival.time','iphase');

        orig_x=zeros(length(orig_lat),1);
        orig_y=zeros(length(orig_lat),1);
        orig_z=zeros(length(orig_lat),1);
        sta_x=zeros(length(orig_lat),1);
        sta_y=zeros(length(orig_lat),1);
        sta_z=zeros(length(orig_lat),1);
        time=zeros(length(orig_lat),1);

        for i=1:length(orig_lat)
                [orig_x(i),orig_y(i),orig_z(i)]=ray_latlon2xyz(orig_lat(i),orig_lon(i),-1*orig_depth(i),ref_lat,ref_lon,ref_alt,projection);
                [sta_x(i),sta_y(i),sta_z(i)]=ray_latlon2xyz(sta_lat(i),sta_lon(i),sta_elev(i),ref_lat,ref_lon,ref_alt,projection);
                time(i)=arr_time(i)-orig_time(i);
                if orig_z(i)>0
                    if strcmp(phase(i),'P')
                        fprintf(fid,'%f %f %f %f %f %f 0 0 %f 1 0\n',orig_x(i),orig_y(i),orig_z(i),sta_x(i),sta_y(i),sta_z(i),time(i));
                    elseif strcmp(phase(i),'S')
                        fprintf(fid,'%f %f %f %f %f %f 0 0 %f 1.73 0\n',orig_x(i),orig_y(i),orig_z(i),sta_x(i),sta_y(i),sta_z(i),time(i));
                    end
                end
        end
    
    case 2
        db = varargin{1};
        locations = varargin{2};
        
        [orid xs ys zs t0 sigx sigy sigz sigt ndat iter rms]=textread(locations,'%d %f %f %f %f %f %f %f %f %d %d %f','headerlines',1);
        
        clear sig* ndat iter rms
        
        db = dbsort(db,'orid');
        
        for i=1:length(orid)
            if zs(i)>0
                db_temp=dbsubset(db,['orid=~/' num2str(orid(i)) '/']);
                [sta_lat,sta_lon,sta_elev,arr_time,phase]=dbgetv(db_temp,'site.lat','site.lon','site.elev','arrival.time','iphase');
                for j=1:length(sta_lat)
                    [sta_x(j),sta_y(j),sta_z(j)]=ray_latlon2xyz(sta_lat(j),sta_lon(j),sta_elev(j),ref_lat,ref_lon,ref_alt,projection);
                    time(j)=arr_time(j)-t0(i);
                    if strcmp(phase(j),'P')
                        fprintf(fid,'%f %f %f %f %f %f 0 0 %f 1 0\n',xs(i),ys(i),zs(i),sta_x(j),sta_y(j),sta_z(j),time(j));
                    elseif strcmp(phase(j),'S')
                        fprintf(fid,'%f %f %f %f %f %f 0 0 %f 1.73 0\n',xs(i),ys(i),zs(i),sta_x(j),sta_y(j),sta_z(j),time(j));
                    end
                end
            else
                disp(['Orid ' num2str(orid(i)) ' has negative depth and will be skipped']) 
            end
        end
end

fclose(fid);